import pandas as pd
import numpy as np
from scipy.stats import norm
import os
import matplotlib.pyplot as plt

# Black-Scholes Function
def black_scholes(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    if option_type == "call":
        return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == "put":
        return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("Invalid option type. Use 'call' or 'put'.")

# Greeks calculation functions
def calculate_greeks(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    delta = norm.cdf(d1) if option_type == "call" else -norm.cdf(-d1)
    gamma = norm.pdf(d1) / (S * sigma * np.sqrt(T))
    vega = S * norm.pdf(d1) * np.sqrt(T) / 100  # Expressed per 1% change in volatility
    theta = - (S * norm.pdf(d1) * sigma) / (2 * np.sqrt(T))
    if option_type == "call":
        theta -= r * K * np.exp(-r * T) * norm.cdf(d2)
    else:
        theta += r * K * np.exp(-r * T) * norm.cdf(-d2)
    theta /= 365  # Per day
    rho = K * T * np.exp(-r * T) * (norm.cdf(d2) if option_type == "call" else -norm.cdf(-d2)) / 100  # Per 1% change in rates

    return delta, gamma, vega, theta, rho

# Load the Excel file
file_path = r'PUT YOUR PERSONAL DIRECTORY'

# Load the second sheet of the Excel file
data = pd.ExcelFile(file_path)
df = data.parse(data.sheet_names[1])  # Load the second sheet

# Data cleaning
columns_mapping = {
    "Unnamed: 0": "Date",
    "ISP IM Equity": "ISP_Price",
    "UCG IM Equity": "UCG_Price",
    "SAN SM Equity": "SAN_Price",
    "BBVA SM Equity": "BBVA_Price",
    "ACA FP Equity": "ACA_Price",
    "BNP FP Equity": "BNP_Price",
    "EUDR1T Curncy": "RiskFreeRate"
}
df = df.rename(columns=columns_mapping)
df = df[~df["Date"].astype(str).str.contains("DATES", na=False)]
df["Date"] = pd.to_datetime(df["Date"], errors='coerce')
df = df.dropna(subset=["Date", "RiskFreeRate"])
df.set_index("Date", inplace=True)

stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
    df_stock = df[[stock, "RiskFreeRate"]].dropna()
    quarterly_strike_prices = df_stock[stock].resample('Q').mean()
    print(f"\nQuarterly Average Prices (Strike Prices) for {stock}:")
    print(quarterly_strike_prices)

    # Calculation of parameters for Black-Scholes
    current_price = df_stock[stock].iloc[-1]
    strike_price = quarterly_strike_prices.iloc[-1]
    time_to_maturity = 0.25  # Quarter in years 
    #We assumed a strike price #given by the quarterly moving average price of the underlying asset
    risk_free_rate = df_stock["RiskFreeRate"].iloc[-1] / 100  # Convert to decimal form
    volatility = df_stock[stock].pct_change().std() * np.sqrt(252)  # Annualized volatility

    # Calculation of option prices
    call_price = black_scholes(current_price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="call")
    put_price = black_scholes(current_price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="put")

    # Output of results
    print(f"Call Option Price for {stock}: {call_price:.2f}")
    print(f"Put Option Price for {stock}: {put_price:.2f}")

    # Creating the price chart and quarterly strike prices
    plt.figure(figsize=(10, 6))

    # Plot the historical stock price (blue line)
    plt.plot(df_stock.index, df_stock[stock], label=f'Price {stock}', color='blue')

    # Plot the Call and Put option prices
    plt.plot(df_stock.index, [black_scholes(price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="call") for price in df_stock[stock]], 
             label=f'Call Price {stock}', color='green', linestyle='--')
    plt.plot(df_stock.index, [black_scholes(price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="put") for price in df_stock[stock]], 
             label=f'Put Price {stock}', color='red', linestyle='--')

    # Plot the strike price (orange dashed line)
    plt.plot(quarterly_strike_prices.index, quarterly_strike_prices, label='Strike Price (Quarterly Average)', color='orange', linestyle='--')

    # Add title, labels, and legend
    plt.title(f'Trend of {stock} Price and Option Prices')
    plt.xlabel('Date')
    plt.ylabel('Price')
    plt.legend()
    plt.grid(True)

    # Show the chart
    plt.show()

Quarterly Average Prices (Strike Prices) for ISP_Price:
Date
2022-03-31    2.399073
2022-06-30    1.944486
2022-09-30    1.756502
2022-12-31    1.997592
2023-03-31    2.380092
2023-06-30    2.366968
2023-09-30    2.462891
2023-12-31    2.542603
Freq: Q-DEC, Name: ISP_Price, dtype: float64
Call Option Price for ISP_Price: 0.24
Put Option Price for ISP_Price: 0.11

Quarterly Average Prices (Strike Prices) for UCG_Price:
Date
2022-03-31    12.527375
2022-06-30     9.711603
2022-09-30     9.712492
2022-12-31    12.398500
2023-03-31    17.019292
2023-06-30    18.972726
2023-09-30    22.032812
2023-12-31    23.888016
Freq: Q-DEC, Name: UCG_Price, dtype: float64
Call Option Price for UCG_Price: 2.48
Put Option Price for UCG_Price: 1.57

Quarterly Average Prices (Strike Prices) for SAN_Price:
Date
2022-03-31    3.116875
2022-06-30    2.881230
2022-09-30    2.508023
2022-12-31    2.682164
2023-03-31    3.342608
2023-06-30    3.261331
2023-09-30    3.515046
2023-12-31    3.672008
Freq: Q-DEC, Name: SAN_Price, dtype: float64
Call Option Price for SAN_Price: 0.32
Put Option Price for SAN_Price: 0.17

Quarterly Average Prices (Strike Prices) for BBVA_Price:
Date
2022-03-31    5.465016
2022-06-30    4.764548
2022-09-30    4.527530
2022-12-31    5.297406
2023-03-31    6.629092
2023-06-30    6.565645
2023-09-30    7.150738
2023-12-31    8.021778
Freq: Q-DEC, Name: BBVA_Price, dtype: float64
Call Option Price for BBVA_Price: 0.66
Put Option Price for BBVA_Price: 0.38

Quarterly Average Prices (Strike Prices) for ACA_Price:
Date
2022-03-31    12.328250
2022-06-30     9.991905
2022-09-30     9.096636
2022-12-31     9.319141
2023-03-31    10.857677
2023-06-30    11.034161
2023-09-30    11.334369
2023-12-31    11.891270
Freq: Q-DEC, Name: ACA_Price, dtype: float64
Call Option Price for ACA_Price: 1.36
Put Option Price for ACA_Price: 0.28

Quarterly Average Prices (Strike Prices) for BNP_Price:
Date
2022-03-31    58.730234
2022-06-30    50.353333
2022-09-30    46.458182
2022-12-31    50.005313
2023-03-31    59.996615
2023-06-30    57.141774
2023-09-30    58.876769
2023-12-31    58.287302
Freq: Q-DEC, Name: BNP_Price, dtype: float64
Call Option Price for BNP_Price: 6.66
Put Option Price for BNP_Price: 1.78

import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline

# Black-Scholes function to calculate option price
def black_scholes(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    if option_type == "call":
        return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == "put":
        return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("Invalid option type. Use 'call' or 'put'.")

# Function to calculate the Greeks
def calculate_greeks(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    delta = norm.cdf(d1) if option_type == "call" else -norm.cdf(-d1)
    gamma = norm.pdf(d1) / (S * sigma * np.sqrt(T))
    vega = S * norm.pdf(d1) * np.sqrt(T) / 100  # Vega per 1% change in volatility
    theta = - (S * norm.pdf(d1) * sigma) / (2 * np.sqrt(T))
    if option_type == "call":
        theta -= r * K * np.exp(-r * T) * norm.cdf(d2)
    else:
        theta += r * K * np.exp(-r * T) * norm.cdf(-d2)
    theta /= 365  # Convert theta to per-day value
    rho = K * T * np.exp(-r * T) * (norm.cdf(d2) if option_type == "call" else -norm.cdf(-d2)) / 100  # Rho per 1% change in rates

    return delta, gamma, vega, theta, rho

# Load the Excel file
file_path = r'PUT YOUR PERSONAL DIRECTORY'

# Load the second sheet from the Excel file
data = pd.ExcelFile(file_path)
df = data.parse(data.sheet_names[1])  # Load second sheet

# Data cleaning and renaming columns
columns_mapping = {
    "Unnamed: 0": "Date",
    "ISP IM Equity": "ISP_Price",
    "UCG IM Equity": "UCG_Price",
    "SAN SM Equity": "SAN_Price",
    "BBVA SM Equity": "BBVA_Price",
    "ACA FP Equity": "ACA_Price",
    "BNP FP Equity": "BNP_Price",
    "EUDR1T Curncy": "RiskFreeRate"
}
df = df.rename(columns=columns_mapping)
df = df[~df["Date"].astype(str).str.contains("DATES", na=False)]
df["Date"] = pd.to_datetime(df["Date"], errors='coerce')
df = df.dropna(subset=["Date", "RiskFreeRate"])
df.set_index("Date", inplace=True)

# List of stocks to analyze
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]
for stock in stocks:
    df_stock = df[[stock, "RiskFreeRate"]].dropna()
    quarterly_strike_prices = df_stock[stock].resample('Q').mean()

    # Compute parameters for Black-Scholes model
    current_price = df_stock[stock].iloc[-1]
    strike_price = quarterly_strike_prices.iloc[-1]
    time_to_maturity = 0.25  # One quarter (3 months) in years
    risk_free_rate = df_stock["RiskFreeRate"].iloc[-1] / 100  # Convert percentage to decimal
    volatility = df_stock[stock].pct_change().std() * np.sqrt(252)  # Annualized volatility

    # Generate a range of underlying asset prices
    prices = np.linspace(current_price * 0.7, current_price * 1.3, 100)  # Range from 70% to 130% of the current price

    # Calculate Greeks for each price in the range
    deltas, gammas, vegas, thetas, rhos = [], [], [], [], []
    for price in prices:
        delta, gamma, vega, theta, rho = calculate_greeks(price, strike_price, time_to_maturity, risk_free_rate, volatility, option_type="call")
        deltas.append(delta)
        gammas.append(gamma)
        vegas.append(vega)
        thetas.append(theta)
        rhos.append(rho)

    # Plot Greeks against stock price
    plt.figure(figsize=(15, 10))

    # Delta vs Price
    plt.subplot(2, 3, 1)
    plt.plot(prices, deltas, label="Delta", color='blue')
    plt.title(f"Delta vs Price for {stock}")
    plt.xlabel('Underlying Price')
    plt.ylabel('Delta')
    plt.grid(True)

    # Gamma vs Price
    plt.subplot(2, 3, 2)
    plt.plot(prices, gammas, label="Gamma", color='green')
    plt.title(f"Gamma vs Price for {stock}")
    plt.xlabel('Underlying Price')
    plt.ylabel('Gamma')
    plt.grid(True)

    # Vega vs Price
    plt.subplot(2, 3, 3)
    plt.plot(prices, vegas, label="Vega", color='purple')
    plt.title(f"Vega vs Price for {stock}")
    plt.xlabel('Underlying Price')
    plt.ylabel('Vega')
    plt.grid(True)

    # Theta vs Price
    plt.subplot(2, 3, 4)
    plt.plot(prices, thetas, label="Theta", color='red')
    plt.title(f"Theta vs Price for {stock}")
    plt.xlabel('Underlying Price')
    plt.ylabel('Theta')
    plt.grid(True)

    # Rho vs Price
    plt.subplot(2, 3, 5)
    plt.plot(prices, rhos, label="Rho", color='orange')
    plt.title(f"Rho vs Price for {stock}")
    plt.xlabel('Underlying Price')
    plt.ylabel('Rho')
    plt.grid(True)

    plt.tight_layout()
    plt.show()

import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    if option_type == "call":
        return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == "put":
        return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("Invalid option type. Use 'call' or 'put'.")

# Assuming 'df' is already cleaned and contains necessary columns
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]

for stock in stocks:
    df_stock = df[[stock, "RiskFreeRate"]].dropna()
    quarterly_strike_prices = df_stock[stock].resample('Q').mean()

    # Extract parameters
    S = df_stock[stock].iloc[-1]  # Current price
    r = df_stock["RiskFreeRate"].iloc[-1] / 100  # Risk-free rate (decimal)
    sigma = df_stock[stock].pct_change().std() * np.sqrt(252)  # Annualized volatility

    # Generate option prices vs. strike price (K)
    K_values = np.linspace(S * 0.8, S * 1.2, 50)
    option_prices_vs_K_call = [black_scholes(S, K, 1, r, sigma, option_type="call") for K in K_values]
    option_prices_vs_K_put = [black_scholes(S, K, 1, r, sigma, option_type="put") for K in K_values]

    # Generate option prices vs. time to maturity (T)
    T_values = np.linspace(0.01, 2, 50)
    option_prices_vs_T_call = [black_scholes(S, S, T, r, sigma, option_type="call") for T in T_values]
    option_prices_vs_T_put = [black_scholes(S, S, T, r, sigma, option_type="put") for T in T_values]

    # Create a figure with subplots for better layout
    fig, ax = plt.subplots(2, 1, figsize=(12, 12))

    # Plot option prices vs. strike price (K)
    ax[0].plot(K_values, option_prices_vs_K_call, label=f'{stock} Call Price vs Strike Price', color='blue', linestyle='-')
    ax[0].plot(K_values, option_prices_vs_K_put, label=f'{stock} Put Price vs Strike Price', color='red', linestyle='--')
    ax[0].set_title(f'{stock} - Option Price vs Strike Price', fontsize=14)
    ax[0].set_xlabel('Strike Price (K)', fontsize=12)
    ax[0].set_ylabel('Option Price', fontsize=12)
    ax[0].grid(True)
    ax[0].legend(loc='best')
    
    # Plot option prices vs. time to maturity (T)
    ax[1].plot(T_values, option_prices_vs_T_call, label=f'{stock} Call Price vs Time to Maturity', color='green', linestyle='-')
    ax[1].plot(T_values, option_prices_vs_T_put, label=f'{stock} Put Price vs Time to Maturity', color='orange', linestyle='--')
    ax[1].set_title(f'{stock} - Option Price vs Time to Maturity', fontsize=14)
    ax[1].set_xlabel('Time to Maturity (T)', fontsize=12)
    ax[1].set_ylabel('Option Price', fontsize=12)
    ax[1].grid(True)
    ax[1].legend(loc='best')

    # Adjust the layout to prevent overlap and display the plots
    plt.tight_layout()
    plt.show()

import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Function to calculate implied volatility
def implied_volatility(option_price, S, K, T, r, option_type="call", tol=1e-5, max_iter=100):
    sigma = 0.2  # Initial guess for volatility
    for i in range(max_iter):
        price = black_scholes(S, K, T, r, sigma, option_type)
        vega = S * norm.pdf((np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))) * np.sqrt(T)
        if vega < 1e-6:  # Avoid division by zero
            return np.nan
        price_diff = option_price - price

        if abs(price_diff) < tol:
            return sigma

        sigma += price_diff / vega
    return np.nan  # Return NaN if no convergence

# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
    if T <= 0 or sigma <= 0:
        return 0.0  # Return 0 for invalid inputs
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    if option_type == "call":
        return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == "put":
        return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("Invalid option type. Use 'call' or 'put'.")

# Load the data from the existing dataframe
# Assuming 'df' contains the cleaned and structured data as in the previous code
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]

for stock in stocks:
    df_stock = df[[stock, "RiskFreeRate"]].dropna()
    quarterly_strike_prices = df_stock[stock].resample('Q').mean()
    print(f"\nCalculating implied volatility surface for {stock}:")

    # Extract parameters
    S = df_stock[stock].iloc[-1]  # Current price
    r = df_stock["RiskFreeRate"].iloc[-1] / 100  # Risk-free rate (decimal)
    
    # Compute historical volatility (annualized)
    daily_returns = df_stock[stock].pct_change().dropna()
    historical_volatility = daily_returns.std() * np.sqrt(252)  # Annualized volatility

    # Debugging step: print key parameters
    print(f"Stock: {stock}, Current Price: {S}, Risk-free Rate: {r:.4f}, Historical Volatility: {historical_volatility:.4f}")
    print("Quarterly Strike Prices:")
    print(quarterly_strike_prices)

    # Grid for strike prices and time to maturity
    K_values = np.linspace(S * 0.8, S * 1.2, 30)
    T_values = np.linspace(0.1, 2, 30)
    K_grid, T_grid = np.meshgrid(K_values, T_values)

    # Compute implied volatilities for the surface
    IV_surface = np.zeros_like(K_grid)
    for i in range(K_grid.shape[0]):
        for j in range(K_grid.shape[1]):
            market_price = black_scholes(S, K_grid[i, j], T_grid[i, j], r, historical_volatility, option_type="call")
            IV_surface[i, j] = implied_volatility(market_price, S, K_grid[i, j], T_grid[i, j], r, option_type="call")

    # Replace NaN values with a small constant for visualization
    IV_surface = np.nan_to_num(IV_surface, nan=0.01)

    # Plot implied volatility surface
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(K_grid, T_grid, IV_surface, cmap='viridis')
    ax.set_title(f"Implied Volatility Surface for {stock}")
    ax.set_xlabel("Strike Price (K)")
    ax.set_ylabel("Time to Maturity (T)")
    ax.set_zlabel("Implied Volatility")
    plt.show()

Calculating implied volatility surface for ISP_Price:
Stock: ISP_Price, Current Price: 2.6435, Risk-free Rate: 0.0398, Historical Volatility: 0.3212
Quarterly Strike Prices:
Date
2022-03-31    2.399073
2022-06-30    1.944486
2022-09-30    1.756502
2022-12-31    1.997592
2023-03-31    2.380092
2023-06-30    2.366968
2023-09-30    2.462891
2023-12-31    2.542603
Freq: Q-DEC, Name: ISP_Price, dtype: float64

Calculating implied volatility surface for UCG_Price:
Stock: UCG_Price, Current Price: 24.565, Risk-free Rate: 0.0398, Historical Volatility: 0.4156
Quarterly Strike Prices:
Date
2022-03-31    12.527375
2022-06-30     9.711603
2022-09-30     9.712492
2022-12-31    12.398500
2023-03-31    17.019292
2023-06-30    18.972726
2023-09-30    22.032812
2023-12-31    23.888016
Freq: Q-DEC, Name: UCG_Price, dtype: float64

Calculating implied volatility surface for SAN_Price:
Stock: SAN_Price, Current Price: 3.7795, Risk-free Rate: 0.0398, Historical Volatility: 0.3215
Quarterly Strike Prices:
Date
2022-03-31    3.116875
2022-06-30    2.881230
2022-09-30    2.508023
2022-12-31    2.682164
2023-03-31    3.342608
2023-06-30    3.261331
2023-09-30    3.515046
2023-12-31    3.672008
Freq: Q-DEC, Name: SAN_Price, dtype: float64

Calculating implied volatility surface for BBVA_Price:
Stock: BBVA_Price, Current Price: 8.226, Risk-free Rate: 0.0398, Historical Volatility: 0.3158
Quarterly Strike Prices:
Date
2022-03-31    5.465016
2022-06-30    4.764548
2022-09-30    4.527530
2022-12-31    5.297406
2023-03-31    6.629092
2023-06-30    6.565645
2023-09-30    7.150738
2023-12-31    8.021778
Freq: Q-DEC, Name: BBVA_Price, dtype: float64

Calculating implied volatility surface for ACA_Price:
Stock: ACA_Price, Current Price: 12.852, Risk-free Rate: 0.0398, Historical Volatility: 0.2803
Quarterly Strike Prices:
Date
2022-03-31    12.328250
2022-06-30     9.991905
2022-09-30     9.096636
2022-12-31     9.319141
2023-03-31    10.857677
2023-06-30    11.034161
2023-09-30    11.334369
2023-12-31    11.891270
Freq: Q-DEC, Name: ACA_Price, dtype: float64

Calculating implied volatility surface for BNP_Price:
Stock: BNP_Price, Current Price: 62.59, Risk-free Rate: 0.0398, Historical Volatility: 0.3107
Quarterly Strike Prices:
Date
2022-03-31    58.730234
2022-06-30    50.353333
2022-09-30    46.458182
2022-12-31    50.005313
2023-03-31    59.996615
2023-06-30    57.141774
2023-09-30    58.876769
2023-12-31    58.287302
Freq: Q-DEC, Name: BNP_Price, dtype: float64

import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    if option_type == "call":
        return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == "put":
        return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("Invalid option type. Use 'call' or 'put'.")

# Function to calculate the local volatility using Dupire's formula (approximation)
def calculate_local_volatility(S, K, T, r, sigma, epsilon=1e-5):
    # Calculate the Black-Scholes price for the option
    C = black_scholes(S, K, T, r, sigma, option_type="call")

    # Calculate the derivative of the option price with respect to time to maturity (T)
    C_T_plus = black_scholes(S, K, T + epsilon, r, sigma, option_type="call")
    C_T_minus = black_scholes(S, K, T - epsilon, r, sigma, option_type="call")
    dC_dT = (C_T_plus - C_T_minus) / (2 * epsilon)  # Numerical approximation

    # Calculate the second derivative of the option price with respect to the stock price (S)
    C_S_plus = black_scholes(S + epsilon, K, T, r, sigma, option_type="call")
    C_S_minus = black_scholes(S - epsilon, K, T, r, sigma, option_type="call")
    d2C_dS2 = (C_S_plus - 2 * C + C_S_minus) / (epsilon ** 2)  # Numerical approximation

    # Local volatility using the Dupire formula approximation
    local_vol = np.sqrt(2 * dC_dT / (S ** 2 * d2C_dS2))
    return local_vol

# Load the data from the existing dataframe
# Assuming 'df' contains the cleaned and structured data as in the previous code
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]

for stock in stocks:
    df_stock = df[[stock, "RiskFreeRate"]].dropna()
    print(f"\nCalculating Local Volatility Surface for {stock}:")

    # Extract parameters
    S_current = df_stock[stock].iloc[-1]  # Current price
    r = df_stock["RiskFreeRate"].iloc[-1] / 100  # Risk-free rate
    sigma = df_stock[stock].pct_change().std() * np.sqrt(252)  # Volatility

    # Debugging step: print key parameters
    print(f"Stock: {stock}, Current Price: {S_current}, Risk-free Rate: {r:.4f}, Volatility: {sigma:.4f}")

    # Grid for stock prices and time to maturity
    K_values = np.linspace(S_current * 0.8, S_current * 1.2, 30)
    T_values = np.linspace(0.1, 2, 30)
    K_grid, T_grid = np.meshgrid(K_values, T_values)

    # Compute local volatilities for the surface
    local_vol_surface = np.zeros_like(K_grid)
    for i in range(K_grid.shape[0]):
        for j in range(K_grid.shape[1]):
            local_vol_surface[i, j] = calculate_local_volatility(S_current, K_grid[i, j], T_grid[i, j], r, sigma)

    # Plot local volatility surface
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(K_grid, T_grid, local_vol_surface, cmap='viridis')
    ax.set_title(f"Local Volatility Surface for {stock}")
    ax.set_xlabel("Strike Price (K)")
    ax.set_ylabel("Time to Maturity (T)")
    ax.set_zlabel("Local Volatility")
    plt.show()

Calculating Local Volatility Surface for ISP_Price:
Stock: ISP_Price, Current Price: 2.6435, Risk-free Rate: 0.0398, Volatility: 0.3212

Calculating Local Volatility Surface for UCG_Price:
Stock: UCG_Price, Current Price: 24.565, Risk-free Rate: 0.0398, Volatility: 0.4156

Calculating Local Volatility Surface for SAN_Price:
Stock: SAN_Price, Current Price: 3.7795, Risk-free Rate: 0.0398, Volatility: 0.3215

Calculating Local Volatility Surface for BBVA_Price:
Stock: BBVA_Price, Current Price: 8.226, Risk-free Rate: 0.0398, Volatility: 0.3158

Calculating Local Volatility Surface for ACA_Price:
Stock: ACA_Price, Current Price: 12.852, Risk-free Rate: 0.0398, Volatility: 0.2803

Calculating Local Volatility Surface for BNP_Price:
Stock: BNP_Price, Current Price: 62.59, Risk-free Rate: 0.0398, Volatility: 0.3107

import pandas as pd
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Black-Scholes pricing function
def black_scholes(S, K, T, r, sigma, option_type="call"):
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)

    if option_type == "call":
        return S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    elif option_type == "put":
        return K * np.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)
    else:
        raise ValueError("Invalid option type. Use 'call' or 'put'.")

# Function to calculate local volatility using Dupire's formula (numerically approximated)
def calculate_local_volatility(S, K, T, r, sigma, epsilon=1e-5):
    # Calculate the Black-Scholes price for the option
    C = black_scholes(S, K, T, r, sigma, option_type="call")

    # Calculate the derivative of the option price with respect to time to maturity (T)
    C_T_plus = black_scholes(S, K, T + epsilon, r, sigma, option_type="call")
    C_T_minus = black_scholes(S, K, T - epsilon, r, sigma, option_type="call")
    dC_dT = (C_T_plus - C_T_minus) / (2 * epsilon)  # Numerical approximation

    # Calculate the second derivative of the option price with respect to the stock price (S)
    C_S_plus = black_scholes(S + epsilon, K, T, r, sigma, option_type="call")
    C_S_minus = black_scholes(S - epsilon, K, T, r, sigma, option_type="call")
    d2C_dS2 = (C_S_plus - 2 * C + C_S_minus) / (epsilon ** 2)  # Numerical approximation

    # Local volatility using the Dupire formula approximation
    local_vol = np.sqrt(2 * dC_dT / (S ** 2 * d2C_dS2))
    return local_vol

# Function to calculate implied volatility using Newton-Raphson method
def implied_volatility(option_price, S, K, T, r, option_type="call", tol=1e-5, max_iter=100):
    sigma = 0.2  # Initial guess for volatility
    for _ in range(max_iter):
        price = black_scholes(S, K, T, r, sigma, option_type)
        vega = S * norm.pdf((np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))) * np.sqrt(T)
        if vega < 1e-6:
            return np.nan  # Avoid division by zero
        price_diff = option_price - price

        if abs(price_diff) < tol:
            return sigma

        sigma += price_diff / vega
    return np.nan  # Return NaN if no convergence

# Load the data from the existing dataframe
# Assuming 'df' contains the cleaned and structured data as in the previous code
stocks = ["ISP_Price", "UCG_Price", "SAN_Price", "BBVA_Price", "ACA_Price", "BNP_Price"]

for stock in stocks:
    df_stock = df[[stock, "RiskFreeRate"]].dropna()
    print(f"\nCalculating Local vs Implied Volatility for {stock}:")

    # Extract parameters
    S_current = df_stock[stock].iloc[-1]  # Current price
    r = df_stock["RiskFreeRate"].iloc[-1] / 100  # Risk-free rate
    sigma = df_stock[stock].pct_change().std() * np.sqrt(252)  # Volatility

    # Debugging step: print key parameters
    print(f"Stock: {stock}, Current Price: {S_current}, Risk-free Rate: {r:.4f}, Volatility: {sigma:.4f}")

    # Grid for stock prices and time to maturity
    K_values = np.linspace(S_current * 0.8, S_current * 1.2, 30)
    T_values = np.linspace(0.1, 2, 30)
    K_grid, T_grid = np.meshgrid(K_values, T_values)

    # Compute local volatilities and implied volatilities for the surface
    local_vol_surface = np.zeros_like(K_grid)
    implied_vol_surface = np.zeros_like(K_grid)
    for i in range(K_grid.shape[0]):
        for j in range(K_grid.shape[1]):
            market_price = black_scholes(S_current, K_grid[i, j], T_grid[i, j], r, sigma, option_type="call")
            implied_vol_surface[i, j] = implied_volatility(market_price, S_current, K_grid[i, j], T_grid[i, j], r, option_type="call")
            local_vol_surface[i, j] = calculate_local_volatility(S_current, K_grid[i, j], T_grid[i, j], r, sigma)

    # Replace NaN values in the implied volatility surface
    implied_vol_surface = np.nan_to_num(implied_vol_surface, nan=0.01)

    # Plot implied volatility surface
    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(121, projection='3d')
    ax.plot_surface(K_grid, T_grid, implied_vol_surface, cmap='viridis')
    ax.set_title(f"Implied Volatility Surface for {stock}")
    ax.set_xlabel("Strike Price (K)")
    ax.set_ylabel("Time to Maturity (T)")
    ax.set_zlabel("Volatility")

    # Plot local volatility surface
    ax2 = fig.add_subplot(122, projection='3d')
    ax2.plot_surface(K_grid, T_grid, local_vol_surface, cmap='plasma')
    ax2.set_title(f"Local Volatility Surface for {stock}")
    ax2.set_xlabel("Strike Price (K)")
    ax2.set_ylabel("Time to Maturity (T)")
    ax2.set_zlabel("Volatility")

    plt.tight_layout()
    plt.show()

Calculating Local vs Implied Volatility for ISP_Price:
Stock: ISP_Price, Current Price: 2.6435, Risk-free Rate: 0.0398, Volatility: 0.3212

Calculating Local vs Implied Volatility for UCG_Price:
Stock: UCG_Price, Current Price: 24.565, Risk-free Rate: 0.0398, Volatility: 0.4156

Calculating Local vs Implied Volatility for SAN_Price:
Stock: SAN_Price, Current Price: 3.7795, Risk-free Rate: 0.0398, Volatility: 0.3215

Calculating Local vs Implied Volatility for BBVA_Price:
Stock: BBVA_Price, Current Price: 8.226, Risk-free Rate: 0.0398, Volatility: 0.3158

Calculating Local vs Implied Volatility for ACA_Price:
Stock: ACA_Price, Current Price: 12.852, Risk-free Rate: 0.0398, Volatility: 0.2803

Calculating Local vs Implied Volatility for BNP_Price:
Stock: BNP_Price, Current Price: 62.59, Risk-free Rate: 0.0398, Volatility: 0.3107